# ifndef _PPFORM_H_
# define _PPFORM_H_

# include"InterpolationStruct.h"
# include"Polynomial.h"

/**
* @brief 定义根据节点和类型，边值条件实现普通多项式样条插值拟合的求解器类
* 
*/
class ppForm
{
private:
    int type = S32;                   // 使用的插值类型，默认为 S32
    std::vector<struct point> points; // 存储节点信息
    std::vector<Polynomial> polyvec;  // 存储分段的多项式插值信息

    // S32 的边值条件
    int boundary = Nature;  // 边值条件，默认 Nature
    double left = 0;        // 存储 x1 处初值条件
    double right = 0;       // 存储 xN 处初值条件

    void ppFormS01();                 // S01 构造函数
    void ppFormS32Complete();         // 三次函数分段插值 S32 Complete 类型
    void ppFormS32Specified_Nature(); // 三次函数分段插值 S32 Specified 和 Nature 类型
    void ppFormS32NotAKnot();         // 三次函数分段插值 S32 Not-a-kont 类型
    void ppFormS32Periodic();         // 三次函数分段插值 S32 Periodic 类型
    
    void TypeForm();        // 根据给定参数进行插值结构构造
    bool judgeXray() const; // 判断输入数据格式是否正确
    
    int search_x(const double & _x) const; // 寻找 x 所在的多项式

public:
    // 初始化构造
    ppForm();
    // 按照节点信息构造
    ppForm(const int & _type, const std::vector<struct point> & _points, 
           const int & _boundary = Nature, const double & _left = 0, const double & _right = 0);
    // 按照节点向量信息构造
    ppForm(const int & _type, const std::vector<double> & x, const std::vector<double> & y,
           const int & _boundary = Nature, const double & _left = 0, const double & _right = 0);
    // 按照自变量值以及给定拟合函数进行构造
    ppForm(const int & _type, const std::vector<double> & x, const Function & func,
           const int & _boundary = Nature);

    double solve(const double & _x) const;       // 根据给定 x 值计算插值多项式计算值
    double diff(const double & _x) const;        // 根据给定 x 值计算插值多项式导数值
    double second_diff(const double & _x) const; // 根据给定 x 值计算插值多项式二阶导数值

    void PrintPoly() const;               // 按照多项式形式打印存储多项式信息
    void PrintPolyCoff() const;           // 按照系数形式打印存储多项式信息
    void PrintPoint() const;              // 打印节点信息
};

// 判断是否有输入错的数据
bool ppForm::judgeXray() const
{
    bool judge = true;
    if (points.size() == 0)
    {
        judge = false;
        return judge;
    }

    for (int i = 0; i < points.size() - 1; i++)
    {
        if (points[i].x >= points[i+1].x)
        {
            judge = false;
            break;
        }
    }
    return judge;
};

// 根据给定类型构造对应的求解器
void ppForm::TypeForm()
{
    // 先对数组进行排序
    std::sort(points.begin(), points.end(), xUpsort);

    // 判断输入点是否横坐标有重复
    if(!judgeXray())
    {
        std::cout << "Values of x are wrong." << std::endl;
        return;
    }

    // 根据类型选择构造
    switch (type)
    {
    case S01: ppFormS01(); break;
    case S32: 
        switch(boundary)
        {
            case Complete: ppFormS32Complete(); break;
            case Specified: case Nature: ppFormS32Specified_Nature(); break;
            case NotAKnot: ppFormS32NotAKnot(); break;
            case Periodic: ppFormS32Periodic(); break;
            default:
                std::cout << "The boundary condition is false." << std::endl;
                break;
        } break;
    default: 
        std::cout << "The type can not be formed." << std::endl;
        break;
    }
};

// 初始化构造
ppForm::ppForm()
{
    struct point p1{0,0}, p2{1,0};
    points.push_back(p1);
    points.push_back(p2);
    TypeForm();
}

// 按照节点信息构造
ppForm::ppForm(const int & _type, const std::vector<struct point> & _points, 
               const int & _boundary, const double & _left, const double & _right) :
    type{_type}, points{_points}, boundary{_boundary}, left{_left}, right{_right}
{
    TypeForm();
};

// 按照向量信息构造
ppForm::ppForm(const int & _type, const std::vector<double> & x, const std::vector<double> & y,
               const int & _boundary, const double & _left, const double & _right) :
    type{_type}, boundary{_boundary}, left{_left}, right{_right}
{
    // 判断点列是否数据对齐
    if (x.size() != y.size())
    {
        std::cout << "Information of points is not completed." << std::endl;
        return;
    }
    for (int i = 0; i < x.size(); i++)
    {
        struct point node{x[i], y[i]};
        points.push_back(node);
    }
    
    TypeForm();
};

// 根据节点与函数进行构造
ppForm::ppForm(const int & _type, const std::vector<double> & x, const Function & func,
               const int & _boundary) :
    type{_type}, boundary{_boundary}
{
    // Compute the boundary condition of Special case
    if(boundary == Complete)
    {
        left = func.diff(x[0]);
        right = func.diff(x[x.size() - 1]);
    }
    else if(boundary == Specified)
    {
        left = func.second_diff(x[0]);
        right = func.second_diff(x[x.size() - 1]);
    }

    for (int i = 0; i < x.size(); i++)
    {
        struct point node{x[i], func(x[i])};
        points.push_back(node);
    }
    
    TypeForm();
};

// 线性分段插值
void ppForm::ppFormS01()
{
    for (int i = 0; i < points.size() - 1; i++)
    {
        double k = (points[i+1].y - points[i].y) / (points[i+1].x - points[i].x);
        double firstcoeff = k;
        double constant = points[i].y - k * points[i].x;

        std::vector<double> coeff{constant, firstcoeff};
        Polynomial poly(coeff);
        polyvec.push_back(poly);
    }
};

// 三次函数分段插值 S32 Complete 类型
void ppForm::ppFormS32Complete()
{
    // Initialize
    int n = points.size(), i;
    std::vector<double> diff01, diff02, mu, lambda;
    for (i = 0; i < n - 1; i++)
    {
        diff01.push_back((points[i + 1].y - points[i].y) / (points[i + 1].x - points[i].x));

        if(i < n - 2)
        {
            mu.push_back((points[i+1].x - points[i].x) / (points[i+2].x - points[i].x));
            lambda.push_back((points[i+2].x - points[i+1].x) / (points[i+2].x - points[i].x));
        }
        
        if(i == 0)
            diff02.push_back((diff01[0] - left) / (points[1].x - points[0].x));
        else
            diff02.push_back((diff01[i] - diff01[i - 1]) / (points[i + 1].x - points[i - 1].x));
    }
    diff02.push_back((right - diff01[n - 2]) / (points[n - 1].x - points[n - 2].x));

    // Calculate s''(xi)
    Eigen::MatrixXd A2 = Eigen::MatrixXd::Zero(n, n);
    Eigen::VectorXd B2(n);
    for (i = 0; i < n; i++)
    {
        if (i == 0)
        {
            A2(0,0) = 2; A2(0,1) = 1; B2(0) = 6 * diff02[0];
        }
        else if (i == n-1)
        {
            A2(n-1, n-2) = 1; A2(n-1, n-1) = 2; B2(n-1) = 6 * diff02[n-1];
        }
        else
        {
            double mul = 1;
            if (mu[i-1] < lambda[i-1])
                mul = 1.0 / mu[i-1];
            else mul = 1.0 / lambda[i-1];
            A2(i, i-1) = mu[i-1] * mul; A2(i, i) = 2 * mul;
            A2(i, i+1) = lambda[i-1] * mul; B2(i) = 6 * diff02[i] * mul;
        }
    }
    
    Eigen::VectorXd mvec2 = A2.colPivHouseholderQr().solve(B2);

    // Compute the vector of polynomials.
    std::vector<double> coe0{-points[0].x, 1};
    Polynomial poly0(coe0);
    double s0 = (mvec2(1) - mvec2(0)) / (points[1].x - points[0].x);
    poly0 = ((s0/6 * poly0 + mvec2(0)/2) * poly0 + left) * poly0 + points[0].y;
    polyvec.push_back(poly0);

    for (i = 1; i < n - 1; i++)
    {
        std::vector<double> coe{-points[i].x, 1};
        Polynomial poly(coe);
        double sp = (mvec2(i+1) - mvec2(i)) / (points[i+1].x - points[i].x);
        double mp = diff01[i] - 1.0 /6 * (mvec2(i+1) + 2*mvec2(i)) * (points[i+1].x - points[i].x);
        poly = ((sp/6 * poly + mvec2(i)/2) * poly + mp) * poly + points[i].y;
        polyvec.push_back(poly);

    }
};

// 三次函数分段插值 S32 Specified 和 Nature 类型
void ppForm::ppFormS32Specified_Nature()
{
    // Initial value
    int n = points.size(), i;
    std::vector<double> diff01, diff02, mu, lambda;
    for (i = 0; i < n - 1; i++)
    {
        diff01.push_back((points[i + 1].y - points[i].y) / (points[i + 1].x - points[i].x));

        if(i < n - 2)
        {
            mu.push_back((points[i+1].x - points[i].x) / (points[i+2].x - points[i].x));
            lambda.push_back((points[i+2].x - points[i+1].x) / (points[i+2].x - points[i].x));
        }
        
        if (i > 0)
        {
            diff02.push_back((diff01[i] - diff01[i-1]) / (points[i+1].x - points[i-1].x));
        }
    }

    // Calculate s''(xi)
    Eigen::MatrixXd A2 = Eigen::MatrixXd::Zero(n-2, n-2);
    Eigen::VectorXd B2(n-2);
    for (i = 0; i < n - 2; i++)
    {
        double mul = 1;
        if(i == 0)
        {
            mul = 1.0 / lambda[0];
            A2(0,0) = 2 * mul; A2(0,1) = 1;
            B2(0) = (6 * diff02[0] - mu[0] * left) * mul;
        }
        else if(i == n-3)
        {
            mul = 1.0 / mu[n-3];
            A2(n-3, n-4) = 1; A2(n-3, n-3) = 2 * mul;
            B2(n-3) = (6 * diff02[n-3] - lambda[n-3] * right) * mul;
        }
        else
        {
            if(mu[i] > lambda[i]) mul = 1.0 / lambda[i];
            else mul = 1.0 / mu[i];

            A2(i, i-1) = mu[i] * mul; A2(i, i) = 2 * mul; A2(i, i+1) = lambda[i] * mul;
            B2(i) = 6 * diff02[i] * mul;
        }
    }
    
    Eigen::VectorXd M2 = A2.colPivHouseholderQr().solve(B2);

    // Compute the vector of polynomials.
    std::vector<double> coe0{-points[0].x, 1};
    Polynomial poly0(coe0);
    double s0 = (M2(0) - left) / (points[1].x - points[0].x);
    double s00 = diff01[0] - 1.0 / 6 * (M2(0) + 2 * left) * (points[1].x - points[0].x);
    poly0 = ((s0/6 * poly0 + left/2) * poly0 + s00) * poly0 + points[0].y;
    polyvec.push_back(poly0);

    for (i = 1; i < n - 2; i++)
    {
        std::vector<double> coe{-points[i].x, 1};
        Polynomial poly(coe);
        double sp = (M2(i) - M2(i-1)) / (points[i+1].x - points[i].x);
        double mp = diff01[i] - 1.0 / 6 * (M2(i) + 2*M2(i-1)) * (points[i+1].x - points[i].x);
        poly = ((sp/6 * poly + M2(i-1)/2) * poly + mp) * poly + points[i].y;
        polyvec.push_back(poly);
    }

    if(n > 2)
    {
        std::vector<double> coen{-points[n - 2].x, 1};
        Polynomial polyn(coen);
        double s0n = (right - M2(n-3)) / (points[n-1].x - points[n-2].x);
        double sn = diff01[n-2] - 1.0 / 6 * (right + 2*M2(n-3)) * (points[n-1].x - points[n-2].x);
        polyn = ((s0n/6 * polyn + M2(n-3)/2) * polyn + sn) * polyn + points[n-2].y;
        polyvec.push_back(polyn);
    }
};

// 三次函数分段插值 S32 Not-a-knot 类型
void ppForm::ppFormS32NotAKnot()
{
    // Initial value
    int n = points.size(), i;
    std::vector<double> diff01, diff02, mu, lambda;
    for (i = 0; i < n - 1; i++)
    {
        diff01.push_back((points[i + 1].y - points[i].y) / (points[i + 1].x - points[i].x));

        if(i < n - 2)
        {
            mu.push_back((points[i+1].x - points[i].x) / (points[i+2].x - points[i].x));
            lambda.push_back((points[i+2].x - points[i+1].x) / (points[i+2].x - points[i].x));
        }
        
        if (i > 0)
        {
            diff02.push_back((diff01[i] - diff01[i-1]) / (points[i+1].x - points[i-1].x));
        }
    }

    // Solving coefficients
    Eigen::MatrixXd A2 = Eigen::MatrixXd::Zero(n, n);
    Eigen::VectorXd B2(n);
    for (i = 0; i < n; i++)
    {
        double mul = 1.0;

        if (i == 0)
        {
            A2(0,0) = points[2].x - points[1].x; A2(0,1) = - points[2].x + points[0].x;
            A2(0,2) = points[1].x - points[0].x; B2(0) = 0;
            if(A2(0,0) < 1) mul = 1.0 / A2(0,0);
            if(A2(0,2) < 1.0 / mul) mul = 1.0 / A2(0,2);
            A2(0,0) *= mul; A2(0,1) *= mul; A2(0,2) *= mul;
        }
        else if (i == n-1)
        {
            A2(n-1, n-3) = points[n-1].x - points[n-2].x;
            A2(n-1, n-2) = - points[n-1].x + points[n-3].x;
            A2(n-1, n-1) = points[n-2].x - points[n-3].x;
            B2(n-1) = 0;
            if(A2(n-1, n-3) < 1) mul = 1.0 / A2(n-1, n-3);
            if(A2(n-1, n-1) < 1.0 / mul) mul = 1.0 / A2(n-1, n-1);
            A2(n-1, n-3) *= mul; A2(n-1, n-2) *= mul; A2(n-1, n-1) *= mul;
        }
        else
        {
            if (mu[i-1] < lambda[i-1])
                mul = 1.0 / mu[i-1];
            else mul = 1.0 / lambda[i-1];
            A2(i, i-1) = mu[i-1] * mul; A2(i, i) = 2 * mul;
            A2(i, i+1) = lambda[i-1] * mul; B2(i) = 6 * diff02[i-1] * mul;
        }
    }
    
    Eigen::VectorXd M2 = A2.colPivHouseholderQr().solve(B2);
    
    // Compute the vector of polynomials.
    for (i = 0; i < n - 1; i++)
    {
        std::vector<double> coe{-points[i].x, 1};
        Polynomial poly(coe);
        double sp = (M2(i+1) - M2(i)) / (points[i+1].x - points[i].x);
        double mp = diff01[i] - 1.0 /6 * (M2(i+1) + 2*M2(i)) * (points[i+1].x - points[i].x);
        poly = ((sp/6 * poly + M2(i)/2) * poly + mp) * poly + points[i].y;
        polyvec.push_back(poly);
    }
};

// 三次函数分段插值 S32 Periodic 类型
void ppForm::ppFormS32Periodic()
{
    // Judge the periodic conditon.
    if(fabs(points[0].y - points[points.size() - 1].y) > 0.000001)
    {
        std::cout << "Cannot satisfy the periodic condition." << std::endl;
        return;
    }
    points[points.size() - 1].y = points[0].y;

    // Calculate
    // Initial value
    int n = points.size(), i;
    std::vector<double> diff01, diff02, mu, lambda;
    for (i = 0; i < n - 1; i++)
    {
        diff01.push_back((points[i + 1].y - points[i].y) / (points[i + 1].x - points[i].x));

        if(i < n - 2)
        {
            mu.push_back((points[i+1].x - points[i].x) / (points[i+2].x - points[i].x));
            lambda.push_back((points[i+2].x - points[i+1].x) / (points[i+2].x - points[i].x));
        }
        
        if (i > 0)
        {
            diff02.push_back((diff01[i] - diff01[i-1]) / (points[i+1].x - points[i-1].x));
        }
    }

    // Solving coefficients
    Eigen::MatrixXd A2 = Eigen::MatrixXd::Zero(n, n);
    Eigen::VectorXd B2(n);
    for (i = 0; i < n; i++)
    {
        double mul = 1.0;

        if (i == 0)
        {
            A2(0,0) = 1; A2(0, n-1) = -1; B2(0) = 0;
        }
        else if (i == n-1)
        {
            A2(n-1, 1) = points[1].x - points[0].x;
            A2(n-1, 0) = 2 * A2(n-1, 1);
            A2(n-1, n-2) = points[n-1].x - points[n-2].x;
            A2(n-1, n-1) = 2 * A2(n-1, n-2);
            B2(n-1) = 6 * (diff01[0] - diff01[n-2]);
        }
        else
        {
            if (mu[i-1] < lambda[i-1])
                mul = 1.0 / mu[i-1];
            else mul = 1.0 / lambda[i-1];
            A2(i, i-1) = mu[i-1] * mul; A2(i, i) = 2 * mul;
            A2(i, i+1) = lambda[i-1] * mul; B2(i) = 6 * diff02[i-1] * mul;
        }
    }
    
    Eigen::VectorXd M2 = A2.colPivHouseholderQr().solve(B2);
    
    // Compute the vector of polynomial.
    for (i = 0; i < n - 1; i++)
    {
        std::vector<double> coe{-points[i].x, 1};
        Polynomial poly(coe);
        double sp = (M2(i+1) - M2(i)) / (points[i+1].x - points[i].x);
        double mp = diff01[i] - 1.0 /6 * (M2(i+1) + 2*M2(i)) * (points[i+1].x - points[i].x);
        poly = ((sp/6 * poly + M2(i)/2) * poly + mp) * poly + points[i].y;
        polyvec.push_back(poly);
    }
};

// 打印节点信息
void ppForm::PrintPoint() const
{
    for (int i = 0; i < points.size(); i++)
    {
        std::cout << points[i].x << "," << points[i].y << " ";
    }
    std::cout << std::endl;
};

// 按照多项式形式打印存储多项式信息
void ppForm::PrintPoly() const
{
    for (int i = 0; i < polyvec.size(); i++)
    {
        polyvec[i].PrintPol();
    }
};

// 按照系数形式打印存储多项式信息
void ppForm::PrintPolyCoff() const
{
    for (int i = 0; i < polyvec.size(); i++)
    {
        polyvec[i].PrintCoff();
    }
};

// 寻找 x 所在的多项式
int ppForm::search_x(const double & _x) const
{
    int le = 0, ri = points.size() - 1, index = 0;
    if(_x >= points[ri].x)
        index = ri - 1;
    else
    {
        while (le < ri)
        {
            int mid = (le + ri) / 2;
            if (_x == points[mid].x)
            {
                index = mid; break;
            }
            else if (_x < points[mid].x) ri = mid;
            else
            {
                if(_x <= points[mid + 1].x)
                {
                    index = mid; break;
                }
                else le = mid + 1;
            }
        }
        if (le == ri) index = le;
    }

    return index;
}


// 根据输入值计算输出的多项式
double ppForm::solve(const double & _x) const
{
    if(_x < points[0].x || _x > points[points.size() - 1].x + 0.00001)
    {
        std::cout << "Out of range." << std::endl;
        return 0;
    }

    return polyvec[search_x(_x)](_x);
}

// 根据给定 x 值计算插值多项式导数值
double ppForm::diff(const double & _x) const
{
    if(_x < points[0].x || _x > points[points.size() - 1].x + 0.00001)
    {
        std::cout << "Out of range." << std::endl;
        return 0;
    }

    return polyvec[search_x(_x)].diff(_x);
};

// 根据给定 x 值计算插值多项式二阶导数值
double ppForm::second_diff(const double & _x) const
{
    if(_x < points[0].x || _x > points[points.size() - 1].x + 0.00001)
    {
        std::cout << "Out of range." << std::endl;
        return 0;
    }

    return polyvec[search_x(_x)].second_diff(_x);
}

# else
// Do nothing.
# endif